home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / rwvector.lha / RWVector2.1 / src / dmath.cc < prev    next >
C/C++ Source or Header  |  1989-08-18  |  3KB  |  141 lines
  1. /*
  2.  *    Double Precision math function definitions
  3.  *
  4.  *    Copyright (C) 1988, 1989.
  5.  *
  6.  *    Dr. Thomas Keffer
  7.  *    Rogue Wave Associates
  8.  *    P.O. Box 85341
  9.  *    Seattle WA 98145-1341
  10.  *
  11.  *    Permission to use, copy, modify, and distribute this
  12.  *    software and its documentation for any purpose and
  13.  *    without fee is hereby granted, provided that the
  14.  *    above copyright notice appear in all copies and that
  15.  *    both that copyright notice and this permission notice
  16.  *    appear in supporting documentation.
  17.  *    
  18.  *    This software is provided "as is" without any
  19.  *    expressed or implied warranty.
  20.  *
  21.  *
  22.  *    @(#)dmath.cc    2.1    8/18/89
  23.  */
  24.  
  25.  
  26. #include "rw/DComplexVec.h"
  27. #include <math.h>
  28.  
  29. DComplexVec
  30. rootsOfOne(int N, unsigned nterms)
  31. {
  32.   /* Use the recursion formulae
  33.    *  sin(a+b) = sin(a)cos(b) + cos(a)sin(b)
  34.    *  cos(a+b) = cos(a)cos(b) - sin(a)sin(b)
  35.    */
  36.  
  37.   if(!N){
  38.     RWnote("::rootsOfOne()", "Zero'th root requested.");
  39.     RWerror(DEFAULT);
  40.   }
  41.  
  42.   unsigned order = abs(N);
  43.  
  44.   double b = 2.0*M_PI/N;
  45.   register double cosb = cos(b);
  46.   register double sinb = sin(b);
  47.  
  48.   register double cosa = 1.0;        // implies a=0
  49.   register double sina = 0.0;
  50.  
  51.   DComplexVec roots(nterms);
  52.   if(nterms){
  53.     roots(0) = DComplex(cosa, sina);
  54.  
  55.     register double temp;
  56.     for(register int i = 1; i<nterms; i++){
  57.       temp = cosa*cosb - sina*sinb;
  58.       sina = sina*cosb + cosa*sinb;
  59.       cosa = temp;
  60.       roots(i) = DComplex(cosa, sina);
  61.     }
  62.   }
  63.  
  64.   return roots;
  65. }
  66.  
  67. /*
  68.  * Expand a complex conjugate even series into its full series.
  69.  */
  70.  
  71. DComplexVec
  72. expandConjugateEven(const DComplexVec& v)
  73. {
  74.   unsigned N = v.length()-1;
  75.   DComplexVec full(2*N);
  76.  
  77.   // Lower half:
  78.   full.slice(0,N+1,1) = v;
  79.   // Upper half:
  80.   full.slice(N+1,N-1,1) = conj(v.slice(N-1,N-1,-1));
  81.  
  82.   return full;
  83. }
  84.  
  85. /*
  86.  * Expand a complex conjugate odd series into its full series.
  87.  */
  88.  
  89. DComplexVec
  90. expandConjugateOdd(const DComplexVec& v)
  91. {
  92.   unsigned N = v.length()-1;
  93.   DComplexVec full(2*N);
  94.  
  95.   // Lower half:
  96.   full.slice(0,N+1,1) = v;
  97.   // Upper half:
  98.   full.slice(N+1,N-1,1) = conj(-v.slice(N-1,N-1,-1));
  99.  
  100.   return full;
  101. }
  102.  
  103. /*
  104.  * Expand an even series into its full series.
  105.  */
  106.  
  107. DoubleVec
  108. expandEven(const DoubleVec& v)
  109. {
  110.   unsigned N = v.length()-1;
  111.   DoubleVec full(2*N);
  112.  
  113.   // Lower half:
  114.   full.slice(0,N+1,1) = v;
  115.   // Upper half:
  116.   full.slice(N+1,N-1,1) = v.slice(N-1,N-1,-1);
  117.  
  118.   return full;
  119. }
  120.  
  121. /*
  122.  * Expand an odd series into its full series.
  123.  */
  124.  
  125. DoubleVec
  126. expandOdd(const DoubleVec& v)
  127. {
  128.   unsigned N = v.length()+1;
  129.   DoubleVec full(2*N);
  130.  
  131.   // Lower half:
  132.   full.slice(1,N-1,1) = v;
  133.   // Upper half:
  134.   full.slice(N+1,N-1,1) = v.slice(N-2,N-1,-1);
  135.  
  136.   full(0) = 0;
  137.   full(N) = 0;
  138.  
  139.   return full;
  140. }
  141.